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Abstract. Accurate measurement of single-trial responses is key to a definitive use of complex 
electromagnetic and hemodynamic measurements in the investigation of brain dynamics. We 
developed the multiple component, Event-Related Potential (mcERP) approach to single-trial 
response estimation io improve our resolution ot dynamic interactions between neuronal 
ensembles located in different layers within a cortical region and/or in different cortical regions. 
The mcERP model asserts that multiple components defined as stereotypic waveforms comprise 
the stimulus-evoked response and that these components may vary in amplitude and latency 
from trial to trial. Maximum a posteriori (MAP) solutions for the model are obtained by 
iterating a set of equal ions derived from the posterior probability. Our first goal was to use the 
mcERP algorithm to analyze interactions (specifically latency and amplitude correlation) 
between responses in different layers within a cortical region. Thus, we evaluated the model by 
applying the algorithm to synthetic data containing two correlated local components and one 
independent far-field component. Three cases were considered: the local components were 
correlated by an interaction in their single-trial amplitudes, by an interaction in their single-trial 
latencies, or by an interaction in both amplitude and latency. We then analyzed the accuracy 
with which the algorithm estimated the component waveshapes and the single-trial parameters as 
a function of the linearity of each of these relationships. Extensions of these analyses to real 
data are discussed as v eil as ongoing work to incorporate more detailed prior information. 


INTRODUCTION 

Electromagnetic and hemodynamic measurements of brain activity are recorded to 
investigate neurological processing of sensory stimuli with the assumption that these 
signals will yield insight into the brain dynamics of sensation. These recorded signals 
represent the superposition of activity from numerous neural sources on the detector 
array. Conventionally, the subject is exposed to repeated presentations or trials of the 
same stimulus, and the recorded signals are averaged across these trials generating the 


average event-related potential (ERP). Although this technique improves signal-to- 
ratio, averaging implicitly assumes that an invariant, stereotypical waveform 
represents the single-trial ERP. This prevents researchers from asking questions about 
the trial-to-trial variability of the evoked response, about dynamic interactions 
between stimulus-eve ked components, and stimulus-induced modulation of ongoing 
EEG rhythms. 

The multiple component Event-Related Potential (mcERP) model was introduced 
as an alternative technique for examining these signals. The model states that multiple 
components generate the signal and that each of these components may vary in 
amplitude and latency from trial to trial. Furthermore, the mcERP model makes no 
assumptions about the interdependency of the different components. Thus, this model 
estimates the different components generated in response to a stimulus presentation, 
quantifies single-trial variability in the evoked potential (Truccolo 2002), and permits 
exploration of dynamic interactions between related components. Anatomical studies 
show connectivity between different cell populations within and between areas, and 
electrophysiological data based on average responses from awake and anesthetized 
animals also support these ideas (Schroeder et al. 1998; Nowak and Bullier 1997, 
Schmolesky et al. 1998). The mcERP model has the potential to show the 
spatiotemporal profile of activation on a single-trial basis and may reveal information 
about parallel and ser ai processing by sensorineural circuitry. 

An interesting basic circuit suggested by the anatomical and electrophysiological 
data is the projection from the granular lamina to the supragranular laminae of VI. 
The granular lamina of VI (lamina IV) represents the primary thalamorecipient layer 
of cortex in the visual system, and the supragranular laminae of VI receive inputs 
from the granular layer before propagating them to the downstream visual areas. The 
anatomical and electrophysiological evidence suggests that there may exist a 
relationship between the amplitudes and latencies of responses simultaneously 
recorded from these different layers. The mcERP algorithm can address this issue 
because it estimates the components from the recorded signals and it also estimates 
their single-trial characteristics. Correlations in the trial-to-trial variability between 
two components would suggest dynamic coupling between those two 
neuroelectrophysiological phenomena. 

Prior to application of any new model to real data however, it is helpful to 
characterize the perfo rmance of the algorithm with a known set of solutions. In this 
paper, the performance of the mcERP algorithm is evaluated on numerous synthetic 
data sets consisting of three physiological components. Each data set contains two 
hypothetical components that display latency coupling, amplitude coupling, or both 
latency and amplitude coupling. The accuracy of the algorithm in estimating the three 
hypothetical components is examined by comparing the estimated parameters of the 
mcERP model with the original parameters utilized in generating the synthetic data. 
Also, the implications of these simulations are discussed with respect to application of 
the algorithm to real cata and with respect to possible future refinements of the model. 


METHODS 


The mcERP Model and Algorithm 


The mcERP model (Knuth et al. 2001, Knuth et al. in preparation) describes 
recorded neuroelectrophysiologic activity as being generated by a set of neural 
ensembles each responding with a stereotypic activation pattern. The stereotypic 
activity is a temporal pattern of electrical activity that we will described as a 
component s„(t), where s is the component’s waveshape across time t and n is an index 
for a particular component. Second, amplitude and latency variability in the evoked 
response occurs from trial to trial (Truccolo 2002), and the mcERP model modifies the 
component waveshape s„{t ) with a trial-specific amplitude scaling factor a„ r and a 
trial-specific latency shift i nr , where n still denotes a particular component and r 
denotes a specific trial. And third, since multiple detectors (referred to as electrodes 
from now on) are often used, a coupling matrix C mn is added to the model to describe 
the relation between a particular electrode m and a particular component n. This is 
important since several detectors may record a particular component differentially. 
These aspects of the mcERP model can be written mathematically as. 
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where n indexes the N components, C m „ denotes the coupling between the m detector 
and the n' h component, a„ r is the amplitude scaling of the n component during the r 
trial, x nr is the latency shift of the n‘ h component during the r‘ h trial, s„(t) is the 
waveshape of the ri h component, and x] mr (t) is the unpredictable signal recorded in the 
m th detector during the r lh trial. 

The most probable set of parameters that satisfy the mcERP model expressed in (1) 
can be calculated from the maximum a posteriori (MAP) solution of the posterior 
probability. The derivation of the equations governing the most probable parameters 
for (1) begins with Bayes’ Theorem, which states that the posterior probability 
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where />(•) is the probability and / is any prior information. The p{data\mode 1,1) in (2) 
is the likelihood of the data given the model, the p(data\I) is called the evidence or the 
prior probability of the data, p(model\I) is called the prior probability of the model, 
and the left-hand side of (2) is the called the posterior probability, which is the 
probability that a given set of model parameters is consistent with the data and the 
information. Bayes' Theorem describes how the prior probability of the model 
p{model\T) is modified by new information such as data. Substituting the parameters 
of the mcERP model nto (2), the posterior probability becomes 
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where the bold terms indicate the entire set of a particular parameter. Since the 
evidence p(x(t)\T) is constant as the model parameters vary, the posterior probability 
can be written 

p(C, s(t), u, t | x(t), /) a p(x(t) | C, s(/j, a, r, /) /?(C,s(/),a,T 1 1) (4) 

with a proportionality constant equal to the inverse of the evidence. Next, the prior 
probability of the model can be factored and expressed as the product of the individual 
probabilities as follows: 
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The p(a\I) and the p(i\I) are assigned uniform distributions and given ranges that are 
physiologically realizable. A Gaussian likelihood is assigned to (5) based on the 
principle of maximum entropy (Sivia, 1996; Jaynes, unpublished) resulting in 
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where a is the expected variance between the predictions and the mean, p(o\I) is the 
prior probability of a, and Q represents the sum of the square of the residual between 
the data and model. {) is derived from (1) as 
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where M is the total number of detectors or electrodes, R is the number of trial 
presentations of a stimulus, T is the number of the data points sampled. A Jeffrey’s 
prior is assigned to a such that p(o\I) = o 1 , and the joint posterior probability of (6) is 
marginalized over all values of a giving the marginal joint posterior probability 

~\mr 

p C, s(/), a, t | x(r), /) oc Q 1 p(C\I)p(s\I). (8) 

The derivation continues with the assignment of prior probability for the coupling 
matrix and the component waveshapes. The solution of the electromagnetic forward 
problem could be used to generate knowledge about the source-detector coupling 
(Knuth, 1998; Knuth & Vaughan, 1999), but for simplicity only uniform prior 
probabilities were assigned to p(C\T) and p(s\I). When is this done and when the 
natural logarithm is taken, the marginal joint posterior probability becomes 

In T = \nQ + const, (9) 

where P is p(C,s(t),a,x\x(t),I). 

The most probable set of parameters for (9) are then determined by solving the 
expression for the maximum a posteriori (MAP) solution of the posterior probability 
by setting the derivative of the posterior with respect to each of the model parameters 
to zero. The resulting equations are not shown here but are found in Knuth et al. 


(2001, in preparation). These equations are then incorporated into an iterative 
algorithm that initially assumes the existence of a single component within the data 
and refines the parameters S;(f), oty r , t ir, 2 nd C m i until there is less than 1% change in 
the waveshape of the s' /(/) or until 15 iterations are complete. After estimating a single 
component, the residual signal is calculated as 

residual m = x mr (t) - £ C mn a nr s n (t -z m ). (1 0) 

n=] 

If the residual displays structure across trials, the mcERP algorithm is applied again to 
estimate a second component S 2 (t) and its related parameters. In this second set of 
iterations, both $ 2 ( 1 ), and their related parameters are updated. This process of 
adding components continues until the user determines that the residual signal 
contains no structure across trials. 

The Synthetic Data 

Synthetic data were created to evaluate the performance of the algorithm on data 
containing two interacting neural sources. The structure of the synthetic data 
simulated field potential recordings from a linear, 15-channel multielectrode array 
spanning laminae 1 through 6 of macaque VI during presentation of a red flash of 
light. In each simulated case, 50 trials of simulated data for each of the 15 channels of 
the electrode array were generated. Three separate components were specified (Figure 
1A): (i) Component 1 resembled activation of a spiny stellate cell population in 



Figure 1. Synthetic data created to evaluate the mcERP algorithm was based on field potential 
recordings from a linear, nultielectrode array positioned to span laminae 1 through 6 simultaneously. 
A. Component I of the synthetic data represents activation of the lamina 4 spiny stellate cells by a 
thalamic input. Compone it 2 represents the subsequent activation of pyramidal cells in the lamina 2 or 
3, and component 3 is a tar-field source whose component is recorded approximately equally by all 
electrode channels. B. The spatial distribution of the three components is observed in their linear 
summation at every electrode channel. Component 1 is located between channels 9 and 10, component 
2 is located between channels 4 and 5, and component 3 is seen nearly equally by each electrode 
channel. 


granular lamina 4 bj a thalamic input; (ii) Component 2 represented activation of 
pyramidal cells in supragranular lamina 2 or 3 by a granular input such as component 
1; and (iii) Component 3 is the electrical activity generated by a far-field source. The 
temporal patterns of these components were chosen to resemble signals observed in 
VI in response to a red flash of light. A coupling matrix specified the relative 
positioning of the different sources such that component 1 was located between 
channels 9 and 10, component 2 was located between channels 4 and 5, and 
component 3 was seen nearly equally across the channels of the electrode array 
(Figure IB). Low-frequency noise with a spectrum of 1/f was added to the data, and 
the standard deviation of the noise was specified as 0.063 so that the signal-to-noise 
ratio (SNR) of the smallest component (component 2) was 1.1 as given by 
SNR=a 2 component 2 nois, where CT denotes the standard deviation across time. 

Components 1 and 2 of the synthetic data were coupled so that component 1 drives 
component 2 via a feedforward connection. We studied four different 
implementations of such a relationship. Method 1 specified a relationship between the 
latency shifts of components 1 and 2 so that when the first component was delayed, so 
was the second, and when the first component was early, so was the second. The 
degree to which the sources were coupled was controlled by a coupling parameter k 
(not to be confused with the coupling matrix). When *=0, the sources are uncoupled 
and the latency of component 2 jitters independent of component 1. When h=l, the 
sources are completely coupled so that component 2 follows component 1 precisely. 
This is written in equation form as 

T 2r =(l-AM0,a,J+*T lr , (11) 


where xj r and X 2 r are the latency shifts of the first and second components respectively 
in the r th trial, £ is a constant that varies the coupling between the components, 
N{0, Gjtat) represents a random number sampled from a normal distribution with mean 
equal to 0 and standard deviation equal to Cjiah and x/r is the latency shift of the first 
component in the / trial. Single-trial synthetic data were generated for the 15 
electrode channels during 50 trials of the stimulus. The amplitude scaling factors of 
the three components in each trial were independent of each other and defined from a 
randomly sampled lognormal distribution with a sample mean \x am p ~ 1 0 and sample 
standard deviation a anp = TO. The latency shifts of component 1 and component 3 
were independent and defined from a normal distribution N{0>10. 0 milliseconds) with 
sample mean \ii a t ~ 0.0 milliseconds and sample standard deviation cj/ a / 10.0 
milliseconds. Eleven separate cases were explored with k = {0.0, 0.1, TO}. In 

each case, Gjtat was equal to 10.0 milliseconds. 

Method 2 specified a coupling between the amplitudes of components 1 and 2. 

This relation was given by: 
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where a/ r and a 2 r are the amplitude scaling factors of the first and second components 
respectively in the t th trial, k is a constant that varies the coupling, N(l,Oj am p ) 
represents a random r umber from a normal distribution with sample mean equal to 0 


and sample standard deviation equal to o jamp , and a j r is the amplitude scaling factor of 
the first component in the r th trial. The first term in the sum on the right-hand side of 
(12) represents the independent lognormal amplitude jitter of component 2, whereas 
the second term represents a sigmoid function whose rise is varied by b (Figure 2). As 
in Method 1, single-trial data for 15 electrode channels were generated across 50 trials. 
The latency shifts of the three components were independent and chosen randomly 
from N(0 t 10.0 milliseconds ). The amplitude scaling factors of components 1 and 3 
were also independen: and chosen from a lognormal distribution with a sample mean 
\x amp = TO and sample standard deviation a amp - TO. Forty different simulations were 
performed using this method with Gj amp ~ 0.75, b - {2°, 2 1 , 2 2 , 2 3 }, and k = {0.1, 0.2, 
1.0} for each value of b. Since a negative value for a is not physiologically 
realizable, each a 2 r value was adjusted by (/-<ct 2 >T where ct 2 represents all amplitude 
scaling factors across R trials and <> represents the mean of that sample. 

Method 3 combined elements of the two previous methods so that the components 
were coupled both by latency shifts and amplitude scaling factors. The latency and 
amplitude relations were specified by (11) and (12) respectively, and the single-trial 
data were again generated for 15 electrode channels during 50 trial presentations of the 
stimulus with the single-trial amplitudes and latencies of components 1 and 3 
determined as before. The latency and amplitude relationships between components 1 
and 2 were calculated with (Sj/at ~ 10.0 milliseconds, ojam P ~ 0.75, b ~ {2 ,2 ,2 ,2 , 
2°, 2 2 , 2 3 }, and k = {0.0, 0.1, ..., TO} for each value of b resulting in 77 simulated 
cases. 

Finally, Method 4 was identical to Method 3 except that the stereotypic waveshape 
of component 2 was altered to be identical to component 1. The latency and 
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Figure 2. Amplitude relation function used in the simulated cases studying amplitude relation alone and 
those studying simultaneous amplitude and latency relation. As the value of b is increased, the relation 
varies from linear to sigmo dal. 



amplitude relations w ere specified with a t i a , - 10.0 milliseconds, a jam p ~ 0.75, b {2 , 
T\ 2°, 2\ 2 2 , 2 3 }, and k = {0.0, 0.1, ..., 1.0} for each value of b. Thus sixty-six 
different combinations were simulated. 

Evaluation of Algorithm Performance 

The mcERP algorithm performance was evaluated with several different measures 
of accuracy. First, the Amari error (Amari et al. 1996) was calculated to quantify the 
algorithm’s ability tc separate each component comprising the synthetic data. The 
basis of this error is derived from the linear mixing model given by 

x = Cs, (13) 

where x is a matrix of the signals recorded in each detector across time, C is the 
mixing (or coupling) matrix describing the component-detector relationship, and s is a 
matrix of the component waveshapes across time. The linear mixing model of (13) 
possesses two indeterm inacies because the coupling coefficients in matrix C and the 
component waveshapes in matrix s can be rescaled by a diagonal scaling matrix E 

x = CIE“'x, (14) 

and the order of the sources can arbitrarily be changed by a permutation matrix n 

x = CSnn _l I _1 5. (15) 

Thus, the linear mixing model can be rewritten as 

x = (CTn)(rr'z-'.s ) (16) 

Expression (16) asserts that estimates of the coupling matrix and the component 
waveshapes will be as accurate as a scaled permutation of the original and may be 
written as 

s = M-'s. (17) 

or 

C = CM. (18) 

where M is a matrix describing the relation between the original matrix and the 
estimated version. A perfect estimate of the coupling matrix or the component 
waveshapes would result in M being a scaled permutation matrix equivalent to A/=EIT. 
Thus, studying M yields insight into the accuracy of the estimated quantity. Solving 
for M l in (17) gives 

M-' ={ss r \ss r y (19) 

and in (18) yields 

AT' =((c T cY(c T cf = (c T c)\c T c), ( 20 ) 

where T denotes the transpose of the marked matrix. The deviation of M 1 from a 
scaled permutation matrix is found by the computing the Amari error 
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where zl we set equal to A f y , i indexes the rows of matrix A,j indexes the columns of 
matrix A, n is the total number of components, k is the index of the position of the 
maximal value, and the denominator normalizes the error. Either (19) or (20) can be 
used to calculate the Amari error; (19) was used to calculate the Amari error for 
relationship methods 1, 2, and 3, and (20) was used to evaluate the Amari error in 
relationship method 4 because the quantity {ss T )'‘ is singular when components 1 and 
2 have the same waveshape. 

Comparing the estimated component waveshapes to the actual component 
waveshapes was the second tool used to evaluate the performance of the algorithm. 
The root-mean-square (RMS) error was calculated as follows: 

r-l - V /=1 ===== 5 (22) 



where j indexes a panicular component, T is the total number of time samples, and s 
denotes the estimated sources s after they have been scaled and permuted so as to be 
comparable to the original sources. 

Accuracy of the latency shifts and amplitude scaling factors was examined by 
calculating the average absolute difference (average ABS error) between the original 
values and the estimated ones. The following two equations define these errors for 
each component j across all R trials: 
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ABS errors below the standard deviation of the variability (Gjamp~l-0 and <5ji at 10.0 
milliseconds) signify an improvement in estimation over that of standard averaging. 


RESULTS 

Performance was first evaluated for synthetic data generated using coupling method 
1. Each measure of performance showed that strong latency coupling does not hamper 
mcERP’s ability to estimate the parameters. For each value of k, the Amari error 
(Figure 3A) remains below 0.06, which indicates less than about 1% error in 



A Latency Relation Alone 

0 1 - • • - 

0 09 ! 
oosl 

I 

0.07 [ 

w 0.06 i 



°0 02 0 4 0 6 0 8 1 

Strength Relation (k) 


B Amplitude Relation Alone 

0 1 - 

; -t- b=l 
: -e- b=2 

0 09 f A b=4 

[ - 0 - b=8 

0 08 
0.07 • 

0.06 


0.05 r 


0 04 1 


0 03 


0 02 " 



Strength of Relation (k) 


Figure 3. The resulting Amari error when the mcERP algorithm is applied to synthetic data with a 
latency relation alone (A) and with an amplitude relation alone (B). 

separation. Although not shown, the component RMS error is below 1 0% and does 
not illustrate a trend with increasing coupling. Finally the average absolute errors of 
the amplitude scaling factors and the latency shifts are well below the standard 
deviation of the original data. 

Method 2 related the amplitude of component 2 to the amplitude of component 1 
via a sigmoid function (12) with a parameter b, which controls the slope of the 
transition (larger b gives a larger slope). Both the transition parameter b (Figure 2) 
and the coupling constant k were varied. Figure 3B illustrates the Amari error for 
method 2. Again the Amari errors were relatively constant as a function of b and k, 
and were less than 0.06. The other parameters were also well estimated. 

Method 3 combined the aspects of methods 1 and 2 such that simultaneous 
amplitude and latency relations were considered. Again, the mcERP algorithm 
estimated the components and their single-trial parameters accurately. Figure 4A 
shows that the Amari error does not vary across k and b. 

The final test of the mcERP algorithm studied synthetic data whose component 1 
and component 2 had the identical waveshapes and were related both in an amplitude 
and latency. Table 1 lists the Amari errors for all of the different cases considered in 
this block of simulations; there are only two cases where the Amari error falls below 
0.06 indicating separation quality was poor in most cases. Figure 4B illustrates a 
surface plot of these values, and shows that as strength of the interaction k increases, 
the Amari error also ncreases. In contrast, as the amplitude interaction varies from 
linear (small b) one sigmoid-shaped (large b), the Amari error does not follow an 
appreciable pattern. Examining the component waveshape, amplitude scaling factor, 
and latency shift errors (all not shown), several different relationships were observed. 

First, component 1 waveshape error increased with increasing k while remaining 
relatively constant wilh varying b. Component l’s estimated amplitude scaling factor 
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Figure 4. mcERP algorith n performance on synthetic data containing both an amplitude and a latency 
relation between components 1 and 2. Panel A shows the Amari error when the related components 
waveshapes differed, and B shows the Amari error when the components' waveshapes were identical. 


error was largest whtn both k and b were small, and estimated latency shift error 
decreased with increasing k while remaining constant with varying b. Second, 
component 2’s RMS waveshape error, average ABS amplitude scaling factor error, 
and average ABS latency shift error all increased with increasing k and all increased 
slightly with increasing b. Finally, the accuracies of component 3 estimates (not 
shown) exhibit no re-ation to both k and b and were equivalent to those in other 
simulations for all practical purposes. 

Samples of the results for two different cases are shown in Figure 5 and Figure 6. 
Figure 5A-C illustra es an example from method 3 where the waveshapes of 
components are different, and where b — 8.0 and k = 1.0. The estimated component 
waveshapes, their spatial distributions, amplitude scaling factors, and latency shifts are 
shown with respect to their original counterparts. Figure 6A-C shows a case from 
relation method 4 where the waveshapes of components 1 and 2 were identical, and 
where b = 8.0 and k = 1.0. In this case, the component waveshapes and their spatial 


Table 1. Amari error when the mcERP algorithm is applied to synthetic data where components 1 and 
2 have the same waveshape and an amplitude and latency relation. 


k value 

! b value 

0.2500 

0.50 

1.00 

2.00 

4.00 

8.00 

0.0 

0.0934 

0.0461 

0.1015 

0.0835 

0.0749 

0.0438 

0.1 

0.0773 

0.0714 

0.0765 

0.0491 

0.0790 

0.0923 

0.2 

0.1081 

0.0771 

0.0942 

0.0888 

0.0966 

0.0482 

0.3 

0.0350 

0.0662 

0.0738 

0.0744 

0.0746 

0.1077 

0.4 

0.0931 

0.0605 

0.1075 

0.1205 

0.1037 

0.1431 

0.5 

0.1160 

0.0949 

0.0930 

0.1407 

0.1379 

0.1452 

0.6 

0.1200 

0.1286 

0.1399 

0.1226 

0.1360 

0.1742 

0.7 

0.1431 

0.1574 

0.1737 

0.1662 

0.1840 

0.2227 

0.8 

0.1521 

0.1664 

0.1848 

0.1805 

0.2104 

0.2241 

0.9 

0.1593 

0.1828 

0.1968 

0.2098 

0.2360 

0.2350 

1.0 

0.1641 

0.1908 

0.2112 

0.2342 

0.2428 

0.2298 



distributions are not well estimated because the components are not separated. In 
addition, method 4 results in greater errors in the estimates of the single-trial features 
as observed in Figure 6B and C in comparison to those illustrated in Figure 5B and C. 


DISCUSSION 

The mcERP model describes the single-trial event-related potential as the linear 
summation of multip e components each described by a stereotypic waveshape and 
possessing a trial-specific amplitude and latency. Application of Bayes’ theorem to 
the model allows development of an algorithm that gives the most probable set of 
parameters to satisfy the model. The advantage of utilizing Bayes’ theorem is that 
algorithm failure results only from model inadequacies, inappropriate or insufficient 
prior information about the situation, or insufficient data. In many cases, model 
inadequacies may be identified by testing the algorithm on synthetic data prior to 
application to real data. 

The performance of the algorithm on data with interacting components is critical as 
dynamical coupling between real cortical components is expected. Furthermore, such 
coupling between components will yield insights into sensory processing within and 
between cortical areas In this paper, the mcERP algorithm was tested using numerous 
synthetic data sets that reconstruct four hypothetical relations between granular and 
supragranular components when evoked by an external stimulus. Amplitude- 
amplitude and latency -latency interactions between two components were found not to 
affect the ability of the algorithm to identify the component waveshapes, their spatial 
locations, or their single-trial amplitude and latency characteristics. These latency- 
latency results were expected as previous studies of the mcERP algorithm found that 
accurate identification of three components with independent amplitude variability but 
no latency variability was possible (Knuth et al. in preparation). This is because a lack 
of latency variability in the components can be described by (1 1) when T/ r = 0 and k = 

1 . 

The second set of simulations explored an amplitude interaction between two 
components. A sigmoid relation was utilized because it provides a threshold and 
saturation. The threshold of the sigmoid relation prevents occurrence of the driven 
component when the driving component’s activation level is below a given threshold. 
The saturation of the sigmoid curve prevents a further increase in the driven 
component’s amplitude when all elements of the neuronal ensemble generating this 
component are “recruited.” In this set of simulations, the algorithm accurately 
estimated all model parameters. 

Simultaneous amplitude and latency interactions between the granular and the 
supragranular components were examined in the third group of simulations. Again, 
the mcERP algorithm accurately estimated all parameters. We hypothesized that even 
though the single-trial characteristics of components 1 and 2 were related, mcERP 
succeeded in identification due to the fact that the waveshapes of the components were 
different. To test this hypothesis, we performed the fourth set of simulations where 
the waveshapes of the coupled components were identical. The quality of separation, 
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Figure 5. mcERP algorithm estimates match the original synthetic data well. In this case, the 
waveshapes of components 1 and 2 are different, k = 1.0 , and b =8.0. A. The component waveshapes 
and their spatial distributions are plotted across electrode channels. The dotted lines in each plot 
indicate the estimated corr ponent, and the solid lines indicate the original waveshapes. The blue lines 
are difficult to observe because the estimated components match the original ones well. B. 
Comparison of the actual and estimated amplitude scaling factors. C. Comparison of the actual and 
estimated latency timeshifs. The dotted diagonal line denotes a perfect match between the actual and 
the estimated values in both B and C. 
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Figure 6. mcERP algorithm estimates do not match the original synthetic data well. In this particular 
case, the waveshapes of components 1 and 2 are equivalent, k — 1.0, and b = 8.0. A. The component 
waveshapes and their spatial distributions are plotted across electrode channels. The dashed lines 
indicate the estimated component, and the solid lines indicate the original waveshapes. B. Comparison 
of the actual and estimated amplitude scaling factors. C. Comparison of the actual and estimated 
latency timeshifts. The dotted, diagonal line denotes a perfect match between the actual and the 
estimated values in both B and C. 


as measured by the Amari error, degrades as the strength of their interaction increases 
but remains relatively unaffected by changes in the shape of the amplitude relation. 
The amplitude scaling factor error and the latency shift error show an inverse 
relationship between the two related components. Since there is considerable mixing 
between these two components, this observation may indicate that the first of the two 
related components estimated by the algorithm is better approximated than the second. 
Indeed, the waveshapes for the highly coupled case in Figure 6 suggest that this is the 
case. Together these findings indicate that a difference in the temporal activation 
pattern of the component is crucial in separating strongly interacting components. 
When the component waveshapes are the same, independent variation of the 
amplitude and latency parameters of the components can assist the algorithm in 
separating the components, but this is not typically sufficient to accurately estimate all 
of the parameters. 

Although several scenarios were considered, other interactions are possible. For 
example, amplitude- atency couplings may exist so that a larger than average 
amplitude of a component in a given trial may result in the earlier activation of a 
second, dynamically coupled component. Investigating the effect of the spatial 
distribution on the estimation of two dynamically coupled sources may comprise a 
second study. Last, the degree to which differences in the waveshapes affect the 
quality of separation might be explored. 

By applying the mcERP algorithm to real data we expect to be able to characterize 
the single-trial properties of these evoked responses, which will provide valuable 
insights into the cortical processing of sensory stimuli. The results of the simulations 
presented here suggest that the mcERP algorithm will be able to distinguish between 
spatially distinct neural ensembles generating differing component waveshapes despite 
the fact that they may be interacting with one another. However, as the waveshapes of 
the two components become more similar, we expect the quality of separation to 
degrade. This scenario may be addressed by incorporating more prior information into 
the mcERP algorithm by adopting a sparse prior on the coupling matrix to enforce 
spatial localization of the sources. Such a prior would also be able to offset the effects 
of real variations in the recording apparatus such as impedance mismatches between 
electrodes. Lastly, we currently employ a crude and rather subjective stopping 
criterion, which could be substantially improved by deriving a criterion that is based 
on the principles of model selection. This would help to avoid possibly both over- 
and under-fitting of :he data. Regardless, the mcERP algorithm is performing as 
expected and the resu ts of these simulations give us confidence the mcERP algorithm 
will produce viable remits when applied to real data. 
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